{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center>\n",
    "<!--img src=\"./img/ods_stickers.jpg\"-->\n",
    "<img src=\"../../img/ods_stickers.jpg\">\n",
    "## Открытый курс по машинному обучению. Сессия № 3\n",
    "\n",
    "### <center> Автор материала: Алексей Борисихин"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <center> Индивидуальный проект по анализу данных </center>\n",
    "### <center> Прогнозирование осадков по данным почасовых наблюдений в Бразилии </center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**План исследования**\n",
    " - Описание набора данных и признаков\n",
    " - Первичный анализ данных\n",
    " - Первичный визуальный анализ данных\n",
    " - Инсайты, найденные зависимости\n",
    " - Выбор метрики\n",
    " - Выбор модели\n",
    " - Предобработка данных\n",
    " - Кросс-валидация и настройка гиперпараметров модели\n",
    " - Создание новых признаков и описание этого процесса\n",
    " - Построение кривых валидации и обучения\n",
    " - Прогноз для тестовой или отложенной выборке\n",
    " - Выводы\n",
    " \n",
    " Более детальное описание <a href=\"https://github.com/Yorko/mlcourse_open/wiki/3rd-session-in-Russian:-all-activities-accounted-for-in-rating-(in-Russian)#%D0%98%D0%BD%D0%B4%D0%B8%D0%B2%D0%B8%D0%B4%D1%83%D0%B0%D0%BB%D1%8C%D0%BD%D1%8B%D0%B9-%D0%BF%D1%80%D0%BE%D0%B5%D0%BA%D1%82\">тут</a>."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-22T19:32:33.174058Z",
     "start_time": "2018-04-22T19:32:31.973149Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.model_selection import train_test_split\n",
    "from scipy import stats\n",
    "\n",
    "import operator\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-22T17:51:34.059022Z",
     "start_time": "2018-04-22T17:51:34.053021Z"
    }
   },
   "outputs": [],
   "source": [
    "# import warnings\n",
    "# warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-22T17:51:34.084071Z",
     "start_time": "2018-04-22T17:51:34.068028Z"
    }
   },
   "outputs": [],
   "source": [
    "ORIGIANL_DATA_FILE = './../data/project_brazil/sudeste.csv'\n",
    "PROCESSED_DATA_FILE = './../data/project_brazil/processed.csv'\n",
    "\n",
    "RANDOM_STATE = 14"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Часть 1. Описание набора данных и признаков"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p>Имеется набор данных почасовых наблюдений погодных условий от 122 погодных станций Юго-Восточного региона Бразилии. В регион входят штаты: Рио-де-Жанейро, Сан-Паулу, Минас-Жерайс и Эспириту-Санту. Данные собираются с 2000 года (однако не все станции начали работы по сбору именно с этого года). Источником данных является INMET (Национальный Метеорологический Институт - Бразилия). Описание и ссылку на скачивание данных можно найти на странице <a href=\"https://www.kaggle.com/PROPPG-PPG/hourly-weather-surface-brazil-southeast-region\">Kaggle Datasets</a>."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Описания признаков в датасете**<br><br>\n",
    "<b>wsid</b> (Numeric) - Уникальный идентификатор погодной станции<br>\n",
    "<b>wsnm</b> (String) - Наименование погодной станции<br>\n",
    "<b>elvt</b> (Numeric) - Высота расположения погодной станции<br>\n",
    "<b>lat</b> (Numeric) - Широта месторасположения погодной станции (град.)<br>\n",
    "<b>lon</b> (Numericw) - Долгота месторасположения погодной станции (град.)<br>\n",
    "<b>inme</b> (String) - INMET-код станции для области<br>\n",
    "<b>city</b> (String) - Город<br>\n",
    "<b>prov</b> (String) - Штат (провинция)<br>\n",
    "<b>mdct</b> (DateTime) - Дата и время наблюдения<br>\n",
    "<b>date</b> (DateTime) - Дата наблюдения<br>\n",
    "<b>yr</b> (Numeric) - Год наблюдения (2000-2016)<br>\n",
    "<b>mo</b> (Numeric) - Месяц наблюдения (1-12)<br>\n",
    "<b>da</b> (Numeric) - День наблюдения (1-31)<br>\n",
    "<b>hr</b> (Numeric) - Час наблюдения (0-23)<br>\n",
    "<b>prcp</b> (Numeric) - Количество осадков за последний час (мм)<br>\n",
    "<b>stp</b> (Numeric) - Давление воздуха (мгновенное) (гПа)<br>\n",
    "<b>smax</b> (Numeric) - Максимальное давление воздуха за последний час (гПа)<br>\n",
    "<b>smin</b> (Numeric) - Минимальное давление воздуха за последний час (гПа)<br>\n",
    "<b>gbrd</b> (String) - Солнечное излучение (кДж/м2)<br>\n",
    "<b>temp</b> (Numeric) - Температура воздуха (мгновенная) (град. Цельсия)<br>\n",
    "<b>dewp</b> (Numeric) - Температура точки росы (мгновенная) (град. Цельсия)<br>\n",
    "<b>tmax</b> (Numeric) - Максимальная температура за последний час (град. Цельсия)<br>\n",
    "<b>dmax</b> (Numericw) - Максимальная температура точки росы за последний час (град. Цельсия)<br>\n",
    "<b>tmin</b> (Numeric) - Минимальная температура за последний час (град. Цельсия)<br>\n",
    "<b>dmin</b> (Numeric) - Минимальная температура точки росыза последний час (град. Цельсия)<br>\n",
    "<b>hmdy</b> (Numeric) - Относительная влажность (мгновенная) (%)<br>\n",
    "<b>hmax</b> (Numeric) - Максимальная относительная влажность за последний час (%)<br>\n",
    "<b>hmin</b> (Numeric) - Минимальная относительная влажность за последний час (%)<br>\n",
    "<b>wdsp</b> (String) - Скорость ветра (м/с)<br>\n",
    "<b>wdct</b> (Numeric) - Направление ветра (град) (0-360)<br>\n",
    "<b>gust</b> (String) - Порывы ветра (м/с)<br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Данный датасет может быть использован для опробования моделей машинного обучения в области прогнозирования погоды. Явно выделенной целевой переменной нет, а на Kaggle предлагается прогнозировать количество осадков или температуру. Мы попробуем спрогнозировать осадки, т.е. целевой меткой будет являться признак **prcp**. Подробнее последует далее.\n",
    "<p>Необходимо учесть тот факт, что данные по осадкам хранятся для того часа, для которого актуальны и прочие наблюдения параметров. Для прогнозирования осадков на час, следующий за наблюдаемыми данными, в тренировочном датасете надо выполнить сдвиг значений целевой переменной на один шаг \"вперед\".\n",
    "<p>Учитывая, что размер датасета значительный (1.72 Гб), сделаем предобработку файла данных, чтобы не загружать весь датасет в память. В простом предварительном взгляде на данные было выявлено, что примерно в 86% исходных данных отсутствует значение признака **prcp**. Выполним фильтрацию таких данных, оставив для дальнейшего анализа и обработки только корректные.\n",
    "<p>Также стоит отметить, что, несмотря на присутствие в данных признаков с датой и временем, датасет не стоит расматривать как упорядоченный (какого-то тренда в погоде с течением времени не может быть в такой короткий временной отрезок, как 10-20 лет). Соответственно, эти признаки несут категориальный характер."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-22T17:52:17.573949Z",
     "start_time": "2018-04-22T17:51:34.089024Z"
    }
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "with open(ORIGIANL_DATA_FILE, 'r', encoding='utf-8') as in_file, open(PROCESSED_DATA_FILE, 'w', encoding='utf-8') as out_file:\n",
    "\n",
    "    total_lines = 0\n",
    "    processed_lines = 0\n",
    "\n",
    "    print('Обработано строк:', end=' ', flush=True)\n",
    "\n",
    "    # read full set\n",
    "    for line in in_file:\n",
    "\n",
    "        total_lines += 1\n",
    "\n",
    "        # write dataset header line\n",
    "        if total_lines == 1:\n",
    "            out_file.write(line)\n",
    "            continue\n",
    "\n",
    "        # progress indication\n",
    "        if total_lines % 1000000 == 0:\n",
    "            print('{}M'.format(total_lines // 1000000), end=' ', flush=True)\n",
    "\n",
    "        # serch 'prcp' field\n",
    "        prcp_field = line.split(',')[14]\n",
    "\n",
    "        if not prcp_field:\n",
    "            continue\n",
    "            \n",
    "        # write line\n",
    "        processed_lines += 1\n",
    "        \n",
    "        out_file.write(line)\n",
    "\n",
    "    print()\n",
    "    print('Всего обработано строк (с заголовком): {}'.format(total_lines))\n",
    "    print('Отфильтрованные строки: {}'.format(processed_lines))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Часть 2. Первичный анализ данных"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Загрузим предварительно обработанные данные."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-22T17:52:27.832664Z",
     "start_time": "2018-04-22T17:52:17.575061Z"
    }
   },
   "outputs": [],
   "source": [
    "df = pd.read_csv(PROCESSED_DATA_FILE)\n",
    "print('Размер предварительно обработанного датасета: {}'.format(df.shape))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-22T17:52:27.906321Z",
     "start_time": "2018-04-22T17:52:27.837663Z"
    }
   },
   "outputs": [],
   "source": [
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Посмотрим на основную информацию о датасете."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-22T17:52:29.571037Z",
     "start_time": "2018-04-22T17:52:27.906321Z"
    }
   },
   "outputs": [],
   "source": [
    "df.info()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Посмотрим на характеристики распределения категориальных признаков. Сразу приведем целочисленные признаки **wsid, yr, mo, da, hr** (которые являются категориальными) к строковому типу для корректной обработки методом describe()."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-22T17:52:42.645667Z",
     "start_time": "2018-04-22T17:52:29.581043Z"
    }
   },
   "outputs": [],
   "source": [
    "int_features = ['wsid', 'yr', 'mo', 'da', 'hr']\n",
    "\n",
    "for feat in int_features:\n",
    "    df[feat] = df[feat].astype('str')\n",
    "    \n",
    "df.describe(exclude=['float64']).T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p>В выборке представлены 4 провинции и 177 городов. Аномальных значений в дате и времени нет. По описанию и практически совпадающим числам уникальных значений признаков **wsid**, **wsnm** и **inme** можно предположить, что данные признаки дублируют друг друга. Проверим это предположение.\n",
    "<p>Идея состоит в следующем. Используя метод *groupby()*, можно сгруппировать значения одного признака по уникальным значениям другого. Затем с помощью метода *describe()* можно посмотреть на число уникальных значений второго признака в группах. Если на всех позициях стоят единицы, то можно утверждать о полном соответствии в значениях двух признаков."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-22T17:52:44.469933Z",
     "start_time": "2018-04-22T17:52:42.649671Z"
    }
   },
   "outputs": [],
   "source": [
    "wsid_groupedby_wsnm_descr = df.groupby('wsnm')['wsid'].describe()\n",
    "wsid_wsnm_non_unique = wsid_groupedby_wsnm_descr[wsid_groupedby_wsnm_descr['unique'] != 1].shape[0]\n",
    "print('Число несоответствий между признаками wsid и wsnm: {}'.format(wsid_wsnm_non_unique))\n",
    "\n",
    "wsid_groupedby_inme_descr = df.groupby('inme')['wsid'].describe()\n",
    "wsid_inme_non_unique = wsid_groupedby_inme_descr[wsid_groupedby_inme_descr['unique'] != 1].shape[0]\n",
    "print('Число несоответствий между признаками wsid и inme: {}'.format(wsid_inme_non_unique))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Рассмотрим несоответствие детальней. Найдем значение **inme** той станции, по которой происходят расхождения, и посмотрим на уникальные значения координат, соответствующие ей."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-22T17:52:44.901534Z",
     "start_time": "2018-04-22T17:52:44.473935Z"
    }
   },
   "outputs": [],
   "source": [
    "non_unique_inme = wsid_groupedby_inme_descr[wsid_groupedby_inme_descr['unique'] != 1].index[0]\n",
    "\n",
    "print(df[df['inme'] == non_unique_inme]['lat'].value_counts())\n",
    "print()\n",
    "print(df[df['inme'] == non_unique_inme]['lon'].value_counts())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Видно, что координаты этой проблемной станции принимают во всем датасете лишь одно значение. Значит, предположение о дублировании друг друга признаками **wsid, wsnm, inme** оказывается верным."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Выделим вещественные признаки и посмотрим на характеристики их распределения."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-22T17:52:47.726378Z",
     "start_time": "2018-04-22T17:52:44.905523Z"
    }
   },
   "outputs": [],
   "source": [
    "df.describe(include=['float64']).T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Можно сделать предварительные выводы по признакам:<br>\n",
    "- **elvt** (высота) - значения характеристик скорее всего не содержат выбросов из-за ошибок измерений, или их число мало; возможно сделать проверку, используя геосервисы.\n",
    "- **lat, lon** (координаты) - присутствуют выбросы в виде нулевых значений координат (все указанные в выборке регионы не паресекаются экватором или нулевым меридианом)\n",
    "- **prcp** (осадки, целевой признак) - хочется верить, что ошибок измерений нет, иначе у нас вопросы к разметке данных; более 50% данных содержат информацию об отсутствии осадков, а 75% - об отсутствии или незначительных осадках\n",
    "- **stp, smax, smin** (показатели давления) - явно просматриваются выбросы из-за ошибок измерений в виде нулевых значений, причем их более 25%\n",
    "- **gbrd** (солнечное излучение) - с этим признакам вообще все плохо. Почти 50% в выборке у этого признака - пропуски. У оставшихся значений более 50% данных содержат нулевые или близкие к нулю выбросы. Выбросы также наблюдаются в виде максимальных значений, в 6 раз превышающих средний показатель солнечного излучения для жарких пустынных мест (ссылки на справочные источники информации будут даны в конце исследования)\n",
    "- **temp, dewp, tmax, dmax, tmin, dmin** (температурные показатели) - значения скорее всего не содержат значительного числа ошибок измерений\n",
    "- **hmdy, hmax, hmin** (показатели влажности) - более 25% данных содержат ошибки измерений в виде нулевых значений\n",
    "- **wdsp, wdct, gust** (показатели ветра) - значения скорее всего не содержат значительного числа ошибок измерений"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Посомтрим на количество пропусков данных для признаков и отобразим их в порядке убывания."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-22T17:52:49.522730Z",
     "start_time": "2018-04-22T17:52:47.729270Z"
    }
   },
   "outputs": [],
   "source": [
    "null_counts = df.isnull().sum()\n",
    "null_counts[null_counts != 0].sort_values(ascending=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Из описания датасета можно увидеть, что признаки **mdct** (дата и время наблюдения) и **date** (дата наблюдения) избыточны, их полностью заменяют следующие признаки: **date, yr, mo, da, hr** (дата, год, месяц, день, час соответственно). Сконвертируем тип признаков **mdct** и **date** во внутренний тип pandas.datetime и проверим, что значения этих признаков соответствуют друг другу."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-22T17:56:34.954023Z",
     "start_time": "2018-04-22T17:52:49.523843Z"
    }
   },
   "outputs": [],
   "source": [
    "df['mdct'] = df['mdct'].apply(pd.to_datetime)\n",
    "df['date'] = df['date'].apply(pd.to_datetime)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-22T18:02:17.916862Z",
     "start_time": "2018-04-22T17:56:34.954023Z"
    }
   },
   "outputs": [],
   "source": [
    "def datetime_accordance_check(row):\n",
    "    return ((row['yr'] == row['mdct'].year) & (row['yr'] == row['date'].year) & \n",
    "            (row['mo'] == row['mdct'].month) & (row['mo'] == row['date'].month) & \n",
    "            (row['da'] == row['mdct'].day) & (row['da'] == row['date'].day) & \n",
    "            (row['hr'] == row['mdct'].hour))\n",
    "\n",
    "df.apply(datetime_accordance_check, axis=1).value_counts()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Видим, что рассматриваемые признаки полностью соответствуют друг другу. Избыточные признаки можно будет удалить из датасета."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Рассмотрим корреляцию вещественных признаков, исключая целевой. Отобразим распределение и топ-20 коррелирующих пар признаков по двум методам: Пирсона и Спирмана (оценивающие силу линейной и монотонной взаимосвязей соответственно)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-22T18:03:29.229459Z",
     "start_time": "2018-04-22T18:02:17.918939Z"
    }
   },
   "outputs": [],
   "source": [
    "numeric_columns = [col for col in df if (df[col].dtype == 'float64')]\n",
    "\n",
    "pearson_corr = df[numeric_columns].corr(method='pearson')\n",
    "spearman_corr = df[numeric_columns].corr(method='spearman')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-22T18:03:29.290425Z",
     "start_time": "2018-04-22T18:03:29.229459Z"
    }
   },
   "outputs": [],
   "source": [
    "no_target_unique_pearson_corr = dict()\n",
    "no_target_unique_spearman_corr = dict()\n",
    "\n",
    "for (row, r_name) in enumerate(pearson_corr.index):\n",
    "    for (col, c_name) in enumerate(pearson_corr.columns[row+1:]):\n",
    "        if('prcp' in [r_name, c_name]):\n",
    "            continue\n",
    "        no_target_unique_pearson_corr['{}-{}'.format(r_name, c_name)] = abs(pearson_corr.iloc[row, col])\n",
    "        \n",
    "for (row, r_name) in enumerate(spearman_corr.index):\n",
    "    for (col, c_name) in enumerate(spearman_corr.columns[row+1:]):\n",
    "        if('prcp' in [r_name, c_name]):\n",
    "            continue\n",
    "        no_target_unique_spearman_corr['{}-{}'.format(r_name, c_name)] = abs(spearman_corr.iloc[row, col])\n",
    "\n",
    "pearson_corr_df = pd.DataFrame(sorted(no_target_unique_pearson_corr.items(), key=operator.itemgetter(1), reverse=True))\n",
    "pearson_corr_df.columns = ['', 'pearson corr']\n",
    "\n",
    "spearman_corr_df = pd.DataFrame(sorted(no_target_unique_spearman_corr.items(), key=operator.itemgetter(1), reverse=True))\n",
    "spearman_corr_df.columns = ['', 'spearman corr']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-22T18:03:29.620434Z",
     "start_time": "2018-04-22T18:03:29.291990Z"
    }
   },
   "outputs": [],
   "source": [
    "corr_df = pd.concat([pearson_corr_df, spearman_corr_df], axis=1)\n",
    "corr_df.hist()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-22T18:03:29.654426Z",
     "start_time": "2018-04-22T18:03:29.625433Z"
    }
   },
   "outputs": [],
   "source": [
    "corr_df = pd.concat([pearson_corr_df.head(20), spearman_corr_df.head(20)], axis=1)\n",
    "corr_df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Видим, что оба метода на топ-результатах дают очень схожие показатели. Можно предположить, что корреляция между этими признаками есть, и взаимосвязь линейна. Но зависимость не объясняется физической природой признаков."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Рассмотри корреляцию с целевым признаком."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-22T18:03:29.699428Z",
     "start_time": "2018-04-22T18:03:29.658431Z"
    }
   },
   "outputs": [],
   "source": [
    "target_pearson_corr = dict(pearson_corr.apply(abs).loc['prcp', :])\n",
    "target_pearson_corr.pop('prcp')\n",
    "\n",
    "target_spearman_corr = dict(spearman_corr.apply(abs).loc['prcp', :])\n",
    "target_spearman_corr.pop('prcp')\n",
    "\n",
    "target_pearson_corr_df = pd.DataFrame(sorted(target_pearson_corr.items(), key=operator.itemgetter(1), reverse=True))\n",
    "target_pearson_corr_df.columns = ['', 'pearson corr']\n",
    "\n",
    "target_spearman_corr_df = pd.DataFrame(sorted(target_spearman_corr.items(), key=operator.itemgetter(1), reverse=True))\n",
    "target_spearman_corr_df.columns = ['', 'spearman corr']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-22T18:03:29.734710Z",
     "start_time": "2018-04-22T18:03:29.702429Z"
    }
   },
   "outputs": [],
   "source": [
    "target_corr_df = pd.concat([target_pearson_corr_df, target_spearman_corr_df], axis=1)\n",
    "target_corr_df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Тут мы уже не наблюдаем схожесть в значениях коэффициентов корреляции. На оценку по методу Пирсона скорее всего повлияла ненормальность распределения, а по методу Спирмана вероятно, что оценки завышены. В любом случае, можно говорить о некоторой монотонной зависимости целевого признака от ряда других, некоторые объясняются физической природой признаков (например, осадки от влажности и точки росы)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Наконец, посмотрим на распределение целевого признака. Построим график в разных масштабах оси ординат."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-22T18:03:31.576655Z",
     "start_time": "2018-04-22T18:03:29.738712Z"
    }
   },
   "outputs": [],
   "source": [
    "_, axes = plt.subplots(nrows=1, ncols=3, figsize=(15, 5))\n",
    "\n",
    "sns.distplot(df['prcp'], ax=axes[0])\n",
    "\n",
    "sns.distplot(df['prcp'], ax=axes[1])\n",
    "axes[1].set_ylim(0, 0.05)\n",
    "\n",
    "sns.distplot(df['prcp'], ax=axes[2])\n",
    "axes[2].set_ylim(0, 0.01)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-22T06:51:40.519284Z",
     "start_time": "2018-04-22T06:51:40.510286Z"
    }
   },
   "source": [
    "Очевидно, что распределение целевого признака далеко от номрального. Тем не менее, проведем статистические тесты на нормальность и скошенность распределения. Для оценки критерия Шапиро-Уилка в модуле scipy.stats существует метод shapiro(). У него есть ограничение на размер выборки - 5000 (в документации указано, что при выборках большего размера точность p-value не гарантируется). Проведем 10 случайных взятий подвыборки заданного размера и найдем среднее значение p-value на этих разбиениях."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-22T18:03:33.175776Z",
     "start_time": "2018-04-22T18:03:31.576655Z"
    }
   },
   "outputs": [],
   "source": [
    "shapiro_p_values = []\n",
    "\n",
    "for split_num in range(10):\n",
    "    np.random.seed(RANDOM_STATE + split_num)\n",
    "    rand_split = np.random.randint(0, df.shape[0], size=5000)\n",
    "    shapiro_p_values.append(stats.shapiro(df.loc[rand_split, 'prcp'])[1])\n",
    "    \n",
    "print('Среднее значение p-value критерия Шапиро-Уилка на случайных разбиениях: {}'.format(np.mean(shapiro_p_values)))\n",
    "print('Значение p-value теста на скошенность распределения: {}'.format(stats.skewtest(df['prcp'])[1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-22T07:55:15.333586Z",
     "start_time": "2018-04-22T07:55:15.134434Z"
    }
   },
   "source": [
    "Полученное среднее p-value критерия Шапиро-Уилка меньше уровня значимости в 0.05, значит распределение не нормальное. Тест на скошенность распределения также не прошел."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Посмотрим внимательней на правый \"хвост\" распределения."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-22T18:03:33.203890Z",
     "start_time": "2018-04-22T18:03:33.180262Z"
    }
   },
   "outputs": [],
   "source": [
    "df[df['prcp'] > 75]['prcp'].value_counts().iloc[-25:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Осадки интенсивностью в районе 100 мм/час - это очень сильные ливни. И хоть такие наблюдения значимы, но их слишком мало в общей массе для какого-то вклада в модель. Выполним фильтрацию нашей целевой переменной по 99% квантилю, чтобы избавиться от таких экстремальных единичных значений."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-22T18:03:34.406788Z",
     "start_time": "2018-04-22T18:03:33.206507Z"
    }
   },
   "outputs": [],
   "source": [
    "prcp_quantile = df['prcp'].quantile(0.99)\n",
    "prcp_filtered = df[df['prcp'] < prcp_quantile]['prcp']\n",
    "\n",
    "plt.figure(figsize=(15, 5))\n",
    "sns.distplot(prcp_filtered)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-22T08:28:59.720651Z",
     "start_time": "2018-04-22T08:28:59.662651Z"
    }
   },
   "source": [
    "<p>После фильтрации по квантилю обнаружилась следующая особенность целевого признака - он принимает дискретные значения (видимо, связанные с ограничениями на измерения в погодных станциях). \n",
    "<p>Учитывая все эти факты, принимаем решение перейти от казалось бы предполагаемой задачи регрессии к задаче многоклассовой классификации. Подчерпнув информацию о принятой классификации интенсивности дождей (точнее, об отсутствии единой таковой), выделим \"на глаз\" следующие классы целевой переменной.\n",
    "- Отсутствие осадков: $prcp = 0.0$\n",
    "- Небольшой дождь: $prcp\\in(0.0; 1.2]$\n",
    "- Умеренный дождь: $prcp\\in(1.2; 2.5]$\n",
    "- Сильный дождь: $prcp > 2.5$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Часть 3. Первичный визуальный анализ данных"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Начнем визуальный анализ с подгруппы категориальных признаков, отражающих временные характеристики: **yr, mo, da, hr**. Построим гистограммы распределений."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-22T18:03:37.355511Z",
     "start_time": "2018-04-22T18:03:34.406788Z"
    }
   },
   "outputs": [],
   "source": [
    "# переведем тип обратно из str\n",
    "int_features = ['wsid', 'yr', 'mo', 'da', 'hr']\n",
    "for feat in int_features:\n",
    "    df[feat] = df[feat].astype('int64')\n",
    "\n",
    "time_features = ['yr', 'mo', 'da', 'hr']\n",
    "\n",
    "_, axes = plt.subplots(nrows=2, ncols=2, figsize=(14, 8))\n",
    "\n",
    "for feat_idx, feat in enumerate(time_features):\n",
    "    ax = axes[feat_idx // 2, feat_idx % 2]\n",
    "    df[feat].hist(ax=ax)\n",
    "    ax.set_title(feat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Видны следующие факты:\n",
    "- количество измерений на погодных станциях резко возросло в 2008 году\n",
    "- в зимние (напомним, что в южном полушарии зима длится с июня по август) месяцы количество измерений сильно меньше, чем в летние"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Построим графики зависимости среднего числа осадков от года и месяца"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-22T18:03:37.967105Z",
     "start_time": "2018-04-22T18:03:37.358867Z"
    }
   },
   "outputs": [],
   "source": [
    "_, axes = plt.subplots(nrows=1, ncols=2, figsize=(14, 4))\n",
    "\n",
    "df[['prcp', 'yr']].groupby('yr').mean().plot(ax=axes[0])\n",
    "df[['prcp', 'mo']].groupby('mo').mean().plot(ax=axes[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Видны следующие факты:\n",
    "- встречаются как засушливые, так и дождливые года\n",
    "- в летние месяцы среднее число осадков меньше, чем в зимние"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Посмотрим на распределение данных по погодным станциям. В первичном анализе данных мы выяснили, что три варианта идентификаторов станций (**wsid, wsnm, inme**) дублируют друг друга. Будем испольовать признак **wsid**, как целочисленный."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-22T18:03:38.477411Z",
     "start_time": "2018-04-22T18:03:37.967105Z"
    }
   },
   "outputs": [],
   "source": [
    "wsid_int = df['wsid'].astype('int64')\n",
    "print('Уникальные занчения признака wsid: {}'.format(np.unique(wsid_int)))\n",
    "\n",
    "wsid_int.hist(bins=len(np.unique(wsid_int)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Количество данных от разных станций в датасете сильно разнится."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Рассмотрим теперь вещественные признаки. Из первичного анализа данных мы знаем, что в ряде этих признаков присутствуют пропуски. Заполним их соответствующими средними значениями."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-22T18:03:38.918104Z",
     "start_time": "2018-04-22T18:03:38.481411Z"
    }
   },
   "outputs": [],
   "source": [
    "real_features = [feat for feat in df if (df[feat].dtype == 'float64')]\n",
    "real_features = list(set(real_features) - {'lat', 'lon', 'prcp'})\n",
    "\n",
    "df_real_features = df[real_features].copy()\n",
    "df_real_features.fillna(df_real_features.mean(), inplace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Построим гистограммы распределений и \"ящики с усами\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-22T18:03:49.622829Z",
     "start_time": "2018-04-22T18:03:38.922100Z"
    }
   },
   "outputs": [],
   "source": [
    "_, axes = plt.subplots(nrows=6, ncols=3, figsize=(16, 20))\n",
    "\n",
    "for feat_idx, feat in enumerate(df_real_features):\n",
    "    ax = axes[feat_idx // 3, feat_idx % 3]\n",
    "    #df_real_features[feat].hist(ax=ax)\n",
    "    sns.distplot(df_real_features[feat], ax=ax)\n",
    "    #ax.set_title(feat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-22T18:03:53.764625Z",
     "start_time": "2018-04-22T18:03:49.622829Z"
    }
   },
   "outputs": [],
   "source": [
    "_, axes = plt.subplots(nrows=6, ncols=3, figsize=(16, 20))\n",
    "\n",
    "for feat_idx, feat in enumerate(df_real_features):\n",
    "    ax = axes[feat_idx // 3, feat_idx % 3]\n",
    "    sns.boxplot(df_real_features[feat], ax=ax)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Видны следующие факты:\n",
    "- как было отмечено в первичном анализе данных, в признаках есть много выбросов с нулевым значением. На гистограммах распределения это видно почти для всех вещественных признаков\n",
    "- без выбросов большая часть признаков близка к нормальному распределению"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Рассмотрим теперь графически корреляцию вещественных признаков между собой и с целевым признаком по методу Спирмана."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-22T18:04:56.525033Z",
     "start_time": "2018-04-22T18:03:53.768626Z"
    }
   },
   "outputs": [],
   "source": [
    "df_real_features_with_target = pd.concat([df_real_features, df['prcp']], axis=1)\n",
    "\n",
    "corr = df_real_features_with_target.corr(method='spearman')\n",
    "\n",
    "# маска для отсечения верхнего треугольника корреляционной матрицы\n",
    "corr_mask = np.zeros_like(corr, dtype=np.bool)\n",
    "corr_mask[np.triu_indices_from(corr_mask)] = True\n",
    "\n",
    "plt.figure(figsize=(10, 10))\n",
    "sns.heatmap(corr, mask=corr_mask, annot=True, fmt='0.1f', linewidths=0.5, cmap='YlGnBu')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Как было отмечено ранее, видна достаточно высокая монотонная корреляция между всеми признаками. Это не очень хорошо, но возможно это связано с большим числом нулевых значений во всех признаках. Также есть признаки, коррелирующие по своей природе. Например, мгновенные, максимальные и минимальные значения какого-то показателя явно будут сильно коррелирующими."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Построим графики попарных зависимостей коррелирующих по своей природе признаков. Это мгновенные, минимальные и максимальные значения для: давления, температуры, температуры точки росы, влажности."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-22T18:06:30.296796Z",
     "start_time": "2018-04-22T18:04:56.525033Z"
    }
   },
   "outputs": [],
   "source": [
    "pairplot_features = ['stp', 'smax', 'smin']\n",
    "\n",
    "sns.pairplot(df_real_features_with_target[pairplot_features])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-22T18:08:03.605872Z",
     "start_time": "2018-04-22T18:06:30.296796Z"
    }
   },
   "outputs": [],
   "source": [
    "pairplot_features = ['temp', 'tmax', 'tmin']\n",
    "\n",
    "sns.pairplot(df_real_features_with_target[pairplot_features])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-22T18:09:37.480127Z",
     "start_time": "2018-04-22T18:08:03.608872Z"
    }
   },
   "outputs": [],
   "source": [
    "pairplot_features = ['dewp', 'dmax', 'dmin']\n",
    "\n",
    "sns.pairplot(df_real_features_with_target[pairplot_features])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-22T18:11:13.780253Z",
     "start_time": "2018-04-22T18:09:37.484899Z"
    }
   },
   "outputs": [],
   "source": [
    "pairplot_features = ['hmdy', 'hmax', 'hmin']\n",
    "\n",
    "sns.pairplot(df_real_features_with_target[pairplot_features])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Исключая из графиков нулевые выбросы, явно просматривается линейный тренд, для всех групп, кроме значений влажности, как это и обнаружилось в первичном анализе данных."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Построим теперь графики попарных зависимостей между целевым признаком и следующими (тут мы исключаем направление ветра и минимальные/максимальные значения температуры, влажности, давления, температуры точки росы): **elvt, stp, gbrd, temp, dewp, hmdy, wdsp, gust**. Для читаемости разделим построение графиков на 2 шага."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-22T18:16:27.080864Z",
     "start_time": "2018-04-22T18:11:13.780253Z"
    }
   },
   "outputs": [],
   "source": [
    "pairplot_features = ['prcp', 'elvt', 'stp', 'gbrd', 'temp']\n",
    "\n",
    "sns.pairplot(df_real_features_with_target[pairplot_features])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-22T18:21:47.660023Z",
     "start_time": "2018-04-22T18:16:27.086600Z"
    }
   },
   "outputs": [],
   "source": [
    "pairplot_features = ['prcp', 'dewp', 'hmdy', 'wdsp', 'gust']\n",
    "\n",
    "sns.pairplot(df_real_features_with_target[pairplot_features])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Явного линейного тренда не видно ни на одной из попарных зависимостей"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Часть 4. Инсайты, найденные зависимости"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Подводя итог первичному анализу данных, сведем воедино все найденные факты и предположения."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Необходимо не забывать о способе получения рассматриваемого датасета - из всех доступных исходных данных были отфильтрованы только те записи, в которых не было пропущено значение целевого признака.\n",
    "- Основная проблема рассматриваемого датасета - большое число выбросов в виде нулевых значений для почти всех признаков. Причем для некоторых признаков это явные выбросы (например, широта с долготой, влажность воздуха), а для других - с первого взгляда это корректные данные (например, температура). При более детальном рассмотрении, идея о корректности данных ставится под сомнение, т.к. слишком высока их доля отосительно других, ненулевых, значений. Скорее всего это связано с особенностью работы погодных станций, которые не измеряют некоторые показатели. \n",
    "- Необходим более детальный анализ, который не укладывается в рамки этой работы. Для заполнения выбросов корректными значениями (и, собственно, принятия решения - выброс это, или действительное значение) можно использовать данные находящихся рядом погодных станций. Другой способ - кластеризовать станции по географическим координатам и высоте, для некорректных данных брать средние значения в сходных условиях из кластера.\n",
    "- Существую две группы в явном виде дублирующих друг друга признаков. Это идентификаторы погодной станции (**wsid, wsnm, inme**) и параметры даты и времени (**mdct, date, yr, mo, da, hr**).\n",
    "- Из даты и времени есть смысл использовать только месяц как категориальный признак. Засушливость года наврятли окажется полезным признаком для прогноза. А день и час распределены почти равномерно, так что польза их тоже очень сомнительна.\n",
    "- Существуют три группы линейно коррелирующих признаков. Это мгновенные, минимальные и максимальные значения для: давления воздуха (**stp, smax, smin**), темепратуры воздуха (**temp, tmax, tmin**), температуры точки росы (**dewp, dmax, dmin**).\n",
    "- Целевой признак (количество осадков) представлен в виде дискретных значений. Сильный скос распределения к нулевым значениям. Учитывая особенности датасета, это могут быть выбросы, как и для других признаков. Но в этом случае наша задача теряет смысл, поэтому остановимся на идее преобладания дней без осадков.\n",
    "- В данной работе будем учитывать идентификатор погодной станции и координаты сами по себе. Опять же, если выходить за рамки этой работы, то более эффективно было бы выделить географические и высотные кластера погодных станций для более очевидного использования моделями факта о географической близости."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Часть 5. Выбор метрики"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Для выбора метрики сначала посмотрим на баланс классов целевой переменной. Выведем долю каждого класса."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-22T18:21:48.695129Z",
     "start_time": "2018-04-22T18:21:47.664016Z"
    }
   },
   "outputs": [],
   "source": [
    "prcp = df['prcp']\n",
    "\n",
    "def prcp_transform(val):\n",
    "    if val == 0.0:\n",
    "        return 0\n",
    "    elif (val > 0.0) & (val <= 1.2):\n",
    "        return 1\n",
    "    elif (val > 1.2) & (val <=2.5):\n",
    "        return 2\n",
    "    elif (val > 2.5):\n",
    "        return 3\n",
    "    \n",
    "prcp = prcp.apply(prcp_transform)\n",
    "\n",
    "prcp.value_counts(normalize=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p>Можно видеть, что классы сильно несбалансированы. Обычной выбором метрики для задачи классификации с сильным дисбалансом классов является метрика **ROC-AUC**. В отличие от более простых метрик (например, доли верных ответов), ROC-AUC учитывает как TPR (True Positive Rate), так и FPR (False Positive Rate). Значит, она не чувствительна к дисбалансу классов. Также эта метрика позволяет дать оценку качества классификации, основываясь на вероятностых предположениях принадлежности к классу, не привязываяь к какому-то конкретному порогу классификации.\n",
    "<p>Т.к. классификация в нашей задаче многоклассовая, то необходимо выбрать способ усреднения метрики. После изучения возможных способов в документации sklearn, выберем метод **'macro'** (итоговая оценка метрики усредняется без учета веса класса, что дает более пессимистичную оценку, чем методы 'micro' и 'weighted')."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Часть 6. Выбор модели"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Для нашей задачи будем рассматривать следующие модели:\n",
    "- **LogisticRegression**. Простая линейная модель, является хорошим бейзлайном в практически любой задачи классификации. Ее линейные коэффициенты также позволят оценить важность того или иного признака в датасете, с помощью L1-регуляризации можно будет избавиться от линейно зависимых признаков. Из минусов - чувствительна к выбросам, что актуально для нашей задачи. Так что не будем ждать от нее хороших результатов.\n",
    "- **XGBoost**. Характер задачи (относительно небольшое число признаков в масштабах выборки, мало категориальных признаков, которые могут \"раздуть\" датасет, большое число выбросов, к которым не чувствительны \"деревянные\" модели) подразумевает, что результат работы градиентного бустинга на деревьях будет одним из лучших среди всех вариантов. Случайный лес, который тоже подходит для описанных условий, оставим без рассмотрения, хотя он тоже мог бы стать хорошим бейзлайном."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Часть 7. Предобработка данных"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Подготовим датасет для дальнейшей работы моделей. Начнем с целевой переменной. Выполним фильтрацию всего датасета по 99% квантилю целевой перемнной. Выполним сдвиг ее значений на одну позицию вперед, чтобы прогнозировать осадки на час, следующий за наблюдаемым, и выкинем из датасета первую строчку, у которой значение признака **prcp** станет пустым."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-22T19:33:36.039109Z",
     "start_time": "2018-04-22T19:33:33.670909Z"
    }
   },
   "outputs": [],
   "source": [
    "prcp_quantile = df['prcp'].quantile(0.99)\n",
    "df = df[df['prcp'] < prcp_quantile]\n",
    "\n",
    "df['prcp'] = df['prcp'].shift(1)\n",
    "df.drop(index=0, inplace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Выполним кодирование значений целевой переменной и выделим ее в отдельный объект."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-22T19:33:38.259513Z",
     "start_time": "2018-04-22T19:33:37.277434Z"
    }
   },
   "outputs": [],
   "source": [
    "y = df['prcp'].apply(prcp_transform).astype('int')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Создадим словарь обработанных признаков, для дальнейшего их объединения в итоговый датасет."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-22T20:11:40.484165Z",
     "start_time": "2018-04-22T20:11:40.420160Z"
    }
   },
   "outputs": [],
   "source": [
    "processed_features = dict()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p>Т.к. одной из моделей у нас является линейная, то категориальные признаки нам необходимо закодировать, а численные - отмасштабировать.\n",
    "<p>Обозначим схему преобразования данных. В нашей задаче у нас отсутствует тестовая выборка, валидацию моделей будем проводить на отложенной части выборки. Но преобразовывать признаки мы будем на всей имеющейся выборке сразу. В общем случае это неправильно, т.к. из валидационной части может \"просочиться\" информация в обучающую часть (например, учтется масштаб всех значений признаков, будут закодированы все значения категориальных признаков, и т.д.). Мы же закроем глаза на это допущение."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Признак **wsid**. Является категориальным, выполним One-Hot-Encoding с выбрасыванием первой колонки (для того, чтобы не возникало зависимости в новых категориальных признаках). Признаки **wsnm, inme** не используем как дублирующие."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-22T20:11:41.943810Z",
     "start_time": "2018-04-22T20:11:41.573276Z"
    }
   },
   "outputs": [],
   "source": [
    "wsid_ohe = pd.get_dummies(df['wsid'], prefix='wsid_ohe', drop_first=True)\n",
    "processed_features['wsid_ohe'] = wsid_ohe"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Признаки **elvt, lat, lon** являются численными, выполним масштабирование."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-22T20:11:42.764963Z",
     "start_time": "2018-04-22T20:11:42.427187Z"
    }
   },
   "outputs": [],
   "source": [
    "stnd_scaler = StandardScaler()\n",
    "\n",
    "elvt_scaled = stnd_scaler.fit_transform(df[['elvt']])\n",
    "processed_features['elvt_scaled'] = pd.DataFrame(elvt_scaled, columns=['elvt_scaled'])\n",
    "\n",
    "lat_scaled = stnd_scaler.fit_transform(df[['lat']])\n",
    "processed_features['lat_scaled'] = pd.DataFrame(lat_scaled, columns=['lat_scaled'])\n",
    "\n",
    "lon_scaled = stnd_scaler.fit_transform(df[['lon']])\n",
    "processed_features['lon_scaled'] = pd.DataFrame(lon_scaled, columns=['lon_scaled'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Признаки **city, prov** не будем использовать в моделях, как малоинформативные (географическую информацию мы уже добавили в виде широты и долготы), но потенциально \"раздувающие\" датасет (они категориальны и придется кодировать)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Как было обозначено ранее, из признаков даты и времени мы будем использовать только **mo** (месяц). Пояснения см. в п. 4. Несмотря на целочисленность признака, он является конечно же категориальным, так что выполним One-Hot-Encoding."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-22T20:11:43.380753Z",
     "start_time": "2018-04-22T20:11:43.228222Z"
    }
   },
   "outputs": [],
   "source": [
    "mo_ohe = pd.get_dummies(df['mo'], prefix='mo_ohe', drop_first=True)\n",
    "processed_features['mo_ohe'] = mo_ohe"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Для **stp, temp, dewp** (давления, температуры воздуха, температуры точки росы) возьмем только мгновенные их значения (максимальные и минимальные значения линейно зависимы, см. п. 4). Воспользуемся масштабированием. Пустые значения заменим средними."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-22T20:11:44.270581Z",
     "start_time": "2018-04-22T20:11:43.881372Z"
    }
   },
   "outputs": [],
   "source": [
    "stp_scaled = stnd_scaler.fit_transform(df[['stp']])\n",
    "processed_features['stp_scaled'] = pd.DataFrame(stp_scaled, columns=['stp_scaled'])\n",
    "\n",
    "temp_scaled = stnd_scaler.fit_transform(df[['temp']].fillna(df[['temp']].mean()))\n",
    "processed_features['temp_scaled'] = pd.DataFrame(temp_scaled, columns=['temp_scaled'])\n",
    "\n",
    "dewp_scaled = stnd_scaler.fit_transform(df[['dewp']].fillna(df[['dewp']].mean()))\n",
    "processed_features['dewp_scaled'] = pd.DataFrame(dewp_scaled, columns=['dewp_scaled'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Для **gbrd** (солнечное излучение) заполним пропуски и выполним масштабирование."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-22T20:11:44.832283Z",
     "start_time": "2018-04-22T20:11:44.674428Z"
    }
   },
   "outputs": [],
   "source": [
    "gbrd_scaled = stnd_scaler.fit_transform(df[['gbrd']].fillna(df[['gbrd']].mean()))\n",
    "processed_features['gbrd_scaled'] = pd.DataFrame(gbrd_scaled, columns=['gbrd_scaled'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Показатели относительной влажности (мгновенная, максимальная и минимальная) используем все, т.к. линейной зависимости в них не наблюдалось, и данные могут быть полезными для прогнозирования осадков. Воспользуемся масштабированием. Пустые значения заменим средними."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-22T20:11:45.997502Z",
     "start_time": "2018-04-22T20:11:45.626492Z"
    }
   },
   "outputs": [],
   "source": [
    "hmdy_scaled = stnd_scaler.fit_transform(df[['hmdy']])\n",
    "processed_features['hmdy_scaled'] = pd.DataFrame(hmdy_scaled, columns=['hmdy_scaled'])\n",
    "\n",
    "hmin_scaled = stnd_scaler.fit_transform(df[['hmin']].fillna(df[['hmin']].mean()))\n",
    "processed_features['hmin_scaled'] = pd.DataFrame(hmin_scaled, columns=['hmin_scaled'])\n",
    "\n",
    "hmax_scaled = stnd_scaler.fit_transform(df[['hmax']].fillna(df[['hmax']].mean()))\n",
    "processed_features['hmax_scaled'] = pd.DataFrame(hmax_scaled, columns=['hmax_scaled'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Показатели ветра **wdsp, gust** (скорость, порывы) отмасштабируем и заполним пропуски. Направление ветра **wdct** разобьем на 10-градусные сектора, 36-ой сектор заменим нулевым, и выполним One-Hot-Encoding."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-22T20:11:48.825999Z",
     "start_time": "2018-04-22T20:11:46.540819Z"
    }
   },
   "outputs": [],
   "source": [
    "wdsp_scaled = stnd_scaler.fit_transform(df[['wdsp']].fillna(df[['wdsp']].mean()))\n",
    "processed_features['wdsp_scaled'] = pd.DataFrame(wdsp_scaled, columns=['wdsp_scaled'])\n",
    "\n",
    "gust_scaled = stnd_scaler.fit_transform(df[['gust']].fillna(df[['gust']].mean()))\n",
    "processed_features['gust_scaled'] = pd.DataFrame(gust_scaled, columns=['gust_scaled'])\n",
    "\n",
    "wdct_processed = df['wdct'].fillna(df[['wdct']].mean())\n",
    "wdct_processed = wdct_processed.apply(lambda x: x // 10).apply(lambda x: x if x != 36.0 else 0.0)\n",
    "wdct_ohe = pd.get_dummies(wdct_processed, prefix='wdct_ohe', drop_first=True)\n",
    "processed_features['wdct_ohe'] = wdct_ohe"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Наконец, соберем все признаки воедино."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-22T20:11:53.757500Z",
     "start_time": "2018-04-22T20:11:51.215013Z"
    }
   },
   "outputs": [],
   "source": [
    "for name_df, part_df in processed_features.items():\n",
    "    processed_features[name_df] = part_df.reset_index(drop=True)\n",
    "\n",
    "X = pd.concat(processed_features.values(), axis=1)\n",
    "print('Размер обработанной выборки: {}'.format(X.shape))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Выделим обучающую и валидационную части выборки. Ранее уже упоминалось, что упорядочивание выборки по времени не имеет смысла, так что воспользуемся случайным разбиением и выделим для валидации 30% всей выборки. Укажем параметр stratify для сохранения баланса классов в новых выборках."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-22T20:11:59.315498Z",
     "start_time": "2018-04-22T20:11:56.118816Z"
    }
   },
   "outputs": [],
   "source": [
    "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, shuffle=True, stratify=y, random_state=RANDOM_STATE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-22T20:12:00.148076Z",
     "start_time": "2018-04-22T20:12:00.134077Z"
    }
   },
   "outputs": [],
   "source": [
    "print('Размер обучающей выборки:', X_train.shape)\n",
    "print('Размер обучающего вектора целевой переменной:', y_train.shape)\n",
    "print('Размер валидационной выборки:', X_test.shape)\n",
    "print('Размер валидационного вектора целевой переменной:', y_test.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Часть 8. Кросс-валидация и настройка гиперпараметров модели"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Часть 9. Создание новых признаков и описание этого процесса"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Часть 10. Построение кривых валидации и обучения"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Часть 11. Прогноз для тестовой или отложенной выборки"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Часть 12. Выводы "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ссылки"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-20T20:19:03.032006Z",
     "start_time": "2018-04-20T20:19:03.020004Z"
    }
   },
   "source": [
    "1. http://meteocenter.net/forecast/comment.htm<br>\n",
    "Классификация осадков по интенсивности была подобрана \"на глаз\", основываясь на классификации из статьи.<br>\n",
    "2. http://www.zodchii.ws/books/info-756.html<br>\n",
    "\"Количество солнечной энергии, поступающей за год на 1 м2 поверхности Земли, изменяется приблизительно от 3000 МДж/м2 на севере до 8000 MДж/м2 в наиболее  жарких пустынных местах\"<br>\n",
    "Что эквивалентно почасовым значениям: 340 кДж/м2 для северных регионов, 910 кДж/м2 для жарких пустынных регионов<br>"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
